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VT) Summary 

o 15 



Since many environmental processes such as heat waves or precipitation are spatial in extent, 
it is likely that a single extreme event affects several locations and the areal modelling of ex- 



i— ' 16 
> 

Uemes is therefore essential if the spatial dependence of extremes has to be appropriately taken 



into account. This paper proposes a framework for conditional simulations of max-stable pro- 
cesses and give closed forms for Brown-Resnick and Schlather processes. We test the method 
on simulated data and give an application to extreme rainfall around Zurich and extreme temper- 
ature in Switzerland. Results show that the proposed framework provides accurate conditional 
simulations and can handle real-sized problems. 

Some key words: Conditional simulation; Markov chain Monte Carlo; Max-stable process; Precipitation; Regular 
conditional distribution; Temperature. 
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c. dombry, f. eyi-mlnko and m. rlbatet 
1. Introduction 

Max-stable processes arise naturally when studying extremes of stochastic 



therefore play a major role in the statistical modelling of spatial extremes (IBuishand et al , 



2008 



Padoan et al 



2010; 



processes and 



Davison et a 



of max-stable proce sses exists dde Haan , 



2011). Although a different spectral characterization 



is (ISchla ther. 



2002) 



1984), for our purposes the most useful representation 



Z(x) = m&xQYi(x), 
i>l 



x € 



(1) 



where {Ci}i>i are the points of a Poisson process on (0, oo) with intensity dA(£) = C 2( K an d 
Yi are independent replicates of a non-negative continuous sample path stochastic process Y 
such that E{Y(x)} = 1 for all x £ M. d . It is well known that Z is a m ax-stable process on R d 



with unit Frechet margins (ide Haan & Fereira , 



2006; 



Schlather, 



20021) . Although <Q} takes the 



pointwise maximum over an infin ite number of po i nts |C? ) and process es Yi, it is possible to get 



approximate realizations from Z (ISchlather! . 



2002; 



Oesting et all 1201 lb - 



Based on CD several 



Brown & Res nick , 



19771 : 



parametric max-stable models have been pr oposed (ISchlather . 



Kabluchko et al. 



2009; 



Davison et al 



2002; 



20111) and share the same finite 



dimensional distribution functions 



pr{Z(xi) < z 1} . . . , Z(x k ) < z k } = exp 



-E 



max 

j=l,...,k 



Y( Xj ) 
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where k € N, Z\ , . . . , z k > an d x-\ x k £ 



Apart from the Smith model (|Genton et al 



functions are explicitly k nown. To bypass this hurdle, 



parametric estimator and IPadoan et al. 
hood estimator. 



20111) , only the bivariate cumul ative distribution 



de Haan & Per eira (2006) propose a semi- 



(|2010r) suggest the use of the maximum pairwise likeli- 
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Conditional simulations of max- stable processes 3 
Paralleling the use of the variogram in class ical geostatistics, the extremal coefficient function 

^3 



( Schlathe r & Tawn , 



2003 



Coolevetal 



2006) 



9(x\ - x 2 ) = -zlogpr{Z(xi) < z,Z(x 2 ) < z] 

is widely used to summarize the spatial dependence of extremes. It takes values in the interval 
[1, 2]; the lower bound indicates complete dependence, and the upper one independence. 
The last decade has see n many advances to develop a geostatistic of extremes and software is 



available to practitioners (IWang . 



2010; 



Schlather, 



2011 



Ribatet, 



20111) . However an important 



tool currently missing is conditional simulation of max-stable processes. In classical geostatis- 



tic ba sed on Gaussian models, conditional simulations are well established (iChiles & Delfinei 



1999) and provide a framework to assess the distribution of a Gaussian random field given val- 



ues observed at fixed locations. For e xample, conditional simulations of Gaussian processes have 



been used to model land topography (IMandelbrot 



1982) 



Condi tiona 



1989, 



1993). 



simulation of max-st able processes is a long-standing problem (IDavis & Resnick 



Wang & Stoevi (12011 



) provide a first solution, but their framework is limited to 
processes having a discrete spectral measure and thus may be too restrictive to appropriately 
model the spatial dependence in complex situations. 

Based on the recent developments on the regular conditional distribution of max-infmitely 
divisible processes, the aim of this paper is to provide a methodology to get conditional sim- 
ulations of max-stable processes with continuous spectral measures. More formally for a study 
region X cM. d , our goal is to derive an algorithm to sample from the regular conditional distribu- 
tion of Z | {Z(xi) = z\, . . . , Z{xk) = Zk} for some z\, . . . , z^ > and k conditioning locations 

xi, . . . ,x k e x. 
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145 2. Conditional simulation of max-stable processes 

146 21. General framework 
147 

This section reviews some key results of an unpublished paper available from the first author 

148 

with a particular emphasis on max-stable processes. Our goal is to give a more practical inter- 
^ pretation of their results from a simulation perspective. To this aim, we recall two key results and 

^ propose a procedure to get conditional realizations of max-stable processes. 

Let be the space of on X C M d and let <P = {<£i}i>i be a Poisson point process 
152 on R x where ^(x) = Q Y i(x) (i = 1,2,.. .) with Ci and Yj as in ®. We write /(x) = 

^ {/(xi), . . . , f(xk)} for all random functions /: X — > M. and x = (xi, . . . , x^) G ^ fc - It is not 

^ difficult to show that for all Borel set ^4 C the Poisson point process {(fi(x)}i>i defined on 

155 jjfc j ias intensity measure 

156 

A X (A) = / pr{CF(x) G^}C 2 dC. 

157 Jo 

158 The point process <P is called regular if the intensity measure A x has an intensity function X x , 

159 i.e., A x (dz) = X x (z)dz, for all x G ^ fc . 

160 The first key point is that provided the point process $ is regular, the intensity function X x and 

161 the conditional intensity function 

162 \(c r \(u, z) , , , 
163 

drives how the conditioning terms {Z(xj) = Zj} (J = 1, . . . , k) are met. 

The second key point is that, conditionally on Z(x) = z, the Poisson point process <3? can be 

165 

decomposed into two independent point processes, say $ = <E>~ U 3> + , where 

166 

16V <j>~ = {tp £ <£> : yj(xj) < zi for alH = 1, . . . , k}, 

^ <P + = {(p G $ : = for some i = 1, . . . , k}. 

169 

170 

171 
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Conditional simulations of max-stable processes 5 
Before introducing a procedure to get conditional realizations of max-stable processes , we 



Wang&Stoevl(l2011 



). 



introduce notation and make connections with the pioneering work of [ 

A function ip G <£ + such that ip(xi) = z^ for some i G {1, . . . , k} is called an extremal func- 
tion associated to xi and denoted by 92+. It is easy to show that there exists almost surely a 
unique extremal function associated to x^. Although <J> + = {92^, . . . , almost surely, it 
might happen that a single extremal function contributes to the random vector Z(x) at sev- 
eral locations xi, e.g., tp£ = tp£ 2 . To take such repetitions into account, we define a random 
partition 9 = (61, . . . , 61) of the set {x±, . . . , x^} into £ = \9\ blocks and extremal functions 
(<pf, ... , tpl) such that <5 + = {iff, . . . ,ip£} and vt(xi) = Zi if G 6j and 92+ (a^) < Z{ if 



0' = 1, k; j = 1, 



Wang & Stoevl (1201 II) call the partition 9 the hitting sce- 



nario. The set of all possible partitions of {x\ , . . . , denoted identifies all possible hitting 
scenarios. 

From a simulation perspective, the fact that <I> _ and <J> + are independent given Z{x) = z 
is especially convenient and suggests a three-step procedure to sample from the conditional 
distribution of Z given Z(x) = z. 

THEOREM 1. Suppose that the point process $ is regular and let (x, s) G X k+m . For r = 
(ri 5 ... ,T£) G 0P k and j = 1,... ,£, define Ij = {i: Xi G tj}, x Tj = (xi) iGlj , z Tj = (zi) i€l] , 
x T c = (xi)i£i. and z T c = (zi)i^j.. Consider the three-step procedure: 

Step 1. Draw a random partition 9 G with distribution 

1 |r| r 

1T X (Z,T) = W {9 = T\ Z(x) = Z} = — -TT X x (Z T .) / K T o\x T .,Z T .(«j) d «i, 

° \ X > Z ) „- =1 J{Uj<Z T c} 3 3 J 

J 3 

where the normalization constant is 

|r| 

C{X,Z)= ^ \[ X > ''-» ) / K f c\ Xf Zfiu^dUj. 
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6 C. Dombry, F. Eyi-Minko and M. Ribatet 

Step 2. Given r = (t\, . . . , ti), draw I independent random vectors ipf(s),..., <pf(s) from 

the distribution 

prj^+fs) G dv | Z(x) = z,6 = rj = ^- ■ J^J l{ u <z T c}\s,x T c)\x Tj ,z Tj (v,u)du^ dv 



where l/.i is the indicator function and 



C 3 = \ 1 {u<z T c}\s,x T c)\x T .,z T .{v,u)dudv, 
J i i 3 3 

and define the random vector Z + (s) = max J=lj £ <f~j(s). 

Step 3. Independently draw a Poisson point process {C«}«>1 on (0, oo) with intensity C -2 dC 
and {Yi}i>\ independent copies of Y, and define the random vector 

Z~(s) = ma^CiYi{s)l{QY,(x)<z}- 

Then the random vector Z(s) = max {Z + (s), Z~ (s)} follows the conditional distribution of 
Z(s) given Z(x) = z. 



The corresponding conditional cumulative distribution function is 

pr{Z(s) < a,Z(x) < z] 



\r\ 



pr{Z(s) < a | Z(x) = z} 



pr{ 



where 



Z(x)<z} E^)II*Ua), (3) 



, . /{ U <s T ^<a} A ( S ,x T c)|^, 2Tj (^^)d^ 

F TJ (a) = pr < a \ Z(x) = z, 6 = r = i — >- — . 

1 3 J J{u<z T c} X x T c\xr r z Tj (u)du 

It is clear from ^ that the conditional random field Z \ {Z(x) = z} is not max-stable. 

2-2. Distribution of the extremal functions 
In this section we derive closed forms for the intensity function \ x (z) and the conditional 
intensity function A a | ^. r ( u) for two widely use d max-stable pro cesses; the Brown-Resnick 



(IBrown & Resnick . 



1977; 



Kabluchko et al 



2009) and the Schlather (Schlather, 



2002) processes. 



Details of the derivations of these closed forms are given in the Appendix. 
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Conditional simulations of max-stable processes 1 
The Brown-Resnick process corresponds to the case where Y(x) = exp{W(x) — j{x)}, x G 

W 1 , in £[]) with W a centered Gaussian process with stationary increments, semi variogram 7 

and such that W(o) = almost surely. For x G X k and provided the covariance matrix of 

the random vector W(x) is positive definite, the intensity function is 

X x (z) = C x exp ^-^ log z T Q x log z + L x log z^j z^ 1 , z G (0, oo) fc , 

with l fc = (l)j=i,...,fe, a\ = {cr 2 (xi)}i=i,...,fc, 



Qx — ^ x 1 



y-ll ,lTy-l 

1^5] x 1^ 



L 3 



1 ( ifS- 1 ^ 



C, = (2 7 r)( 1 - fc )/2| Sx |-i/2 (1 T S -i U , ) -i/2 exp 



^x V^ ^T v-1 
f 1 (l^E-V^ - l) 2 1 ,T„_i 2 



°x ^X CT X f > 



v 2 1^1, 2 
and for all (s, x) G x m+k , (u, z) G (0, oo) m+A: and provided the covariance matrix E/ s x -\ is pos- 
itive definite, the conditional intensity function corresponds to a multivariate log-normal proba- 
bility density function 

f 1 1 m 

K\x,z{ u ) = ( 27r )~ m/2 | s s |xr 1/2 exp I _-(logu- // s | :Ci2 ) T S^(logu - n s \ XtZ ) \ Y[ur\ 

where /i s | Xj2 G M m and E a u are the mean and covariance matrix of the underlying normal dis- 
tribution and are given by 



^s|x ~~ Jm,kQ(s,x)Jm,k, 



Ms|x,z — ^I J {s,x)Jm,k — log Z 1 ' Jm,kQ(s,x)Jm,k^ E 



s|x) 



with 



J, 



m,k 



Id, 



J, 



m,k 



0m,fc 

Id fe 



where Id^ denotes the x k identity matrix and mi k the m x n null matrix. 

The Schlather process considers the case where Y(x) = (2-k) 1 / 2 max{0, e(x)}, x G R d , in © 
with e a standard Gaussian process with correlation function p. The associated point process 
is not regular and it is more convenient to consider the equivalent representation where Y(x) = 
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337 (27r) 1//2 e(x), x £ For x £ X k and provided the covariance matrix H x of the random vector 

338 e(x) is positive definite, the intensity function is 
339 



340 



345 
346 



360 
361 
362 
363 



\ x ( Z ) = n-( k - i y^ x r 1/2 aAzr ik+1]/2 r^ 



k ' 1 ^ z £ R\ 



where a x (z) = 1 z. 

For (s, x) £ X m+k , (u, z) £ M m+fc and provided that the covariance matrix S( s x ) is positive 
definite, the conditional intensity function \ 8 \ x z (u) corresponds to the density of a multivariate 
^44 Student distribution with k + 1 degrees of freedom, location parameter \i = T ls:x T,~ 1 z, and scale 

matrix 



347 S - (S< - Ss^S^ 1 S x:s ) , S(a,x) 



348 
349 
350 
351 
352 
353 
354 
355 
356 
357 
358 

359 pr(0 £ • | Q-j = T-j) t (4) 



3. Markov chain Monte Carlo sampler 

The previous section introduced a procedure to get realizations from the regular conditional 
distribution of max-stable processes. This sampling scheme amounts to sample from a discrete 
distribution whose state space corresponds to all possible partitions of the set of conditioning 
points, see Theorem Q] step 1. Hence, even for a moderate number k of conditioning locations, 
the state space becomes very large and the distribution tt x (z, •) cannot be computed. It turns 
out that a Gibbs sampler is especially convenient. 

For r £ let T-j be the restriction of r to the set {x\, . . . , x^} \ {%j}- Our goal is to 
simulate from the conditional distribution 



where 6 £ 0?^, is a random partition which follows the target distribution ir x (z, •). 



Conditional simulations of max-stable processes 9 

385 Since the number of possible updates is always less than k, a combinatorial explosion is 

386 avoided. Indeed for r 6 of size £, the number of partitions r* G such that r*^ = T_j 

387 for some j e {1, . . . , k} is 



388 
389 
390 
391 
392 
393 
394 
395 
396 



b + 



£ if {xj} is a partitioning set of r, 

£ + 1 if {xj} is not a partitioning set of r, 

since the point x^ may be reallocated to any partitioning set of r_j or to a new one. 

To illustrate consider the set {xi,x 2 ,x 3 } and let r = ({xi, X2}, {^3})- Then the possible 
partitions r* such that r* 2 = r_ 2 are ({xi, x 2 }, {x 3 }), ({xi}, {x 2 }, {x 3 }), ({xi}, {x 2 , x 3 }), 
while there exists only two partitions such that r* 3 = r_ 3 , i.e., ({xi, x 2 }, {x 3 }), ({xi, x 2 , x 3 }). 

The distribution (0]) has nice properties. Since for all r* G such that r* ■ = T_j we have 



pr [^ = r* 1 ^ = r_,] = ■ "\ & 



W T J = X x (Zrj) / X x \ x z {u)du. 

J{u<z T c] i J J 
j 

Since many factors cancel out on the right hand side of ©, the Gibbs sampler is especially 
convenient. 

The most computationally demanding part of ((5]) is the evaluation of the integral 



397 fe^ fe 
29 § where 

399 
400 
401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
411 



/ K T c\x T .,Z T {u)dU. 

J\u<z,c\ i 3 3 



'{u<Z T c} j 
3 

For the Brown-Resnick and Schlather processes, we follow the lines of lGenzl (I1992T) and compute 
these probabilities using a separation of variables method which provides a transformation of the 
original integration problem to the unit hyper-cube. Further a quasi Monte Carlo scheme and 
antithetic variable sampling is used to improve efficiency. 
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Since it is not obvious how to implement a Gibbs sampler whose target distribution has support 

the remainder of this section gives practical details. For any fixed locations x\, . . . , S 

X, we first describe how each partition of {x\, . . . , Xk} is stored. To illustrate consider the set 

{x±, X2, X3} and the partition ({x\, x 2 }; {^3})- This partition is defined as (1, 1, 2), indicating 

that x\ and X2 belong to the same partitioning set labeled 1 and X3 belongs to the partitioning set 

2. There exist several equivalent notations for this partition: for example one can use (2, 2, 1) or 

(1,1,3). Since there is a one-one mapping between 3 d ^ and the set 

= \ • • • > a k) '■ i G {2, . . • , k}, 1 = a\ < a,i < max a,- + 1, a% € Z > , 
[ i<i<« J 

we shall restrict our attention to the partitions that live in 3?^ and going back to our example we 
see that (1, 1, 2) is valid but (2, 2, 1) and (1, 1, 3) are not. 

For r 6 3?^ of size I, let r\ = Yli=i ^{n=aj} an d r 2 = St=i ^{n=b}^ ie > tne number of 
conditioning locations that belong to the partitioning sets aj and b where b € {1, . . . , b + } with 



b+ 



(ri = 1), 
+ 1 (ri^l). 



Then the conditional probability distribution ([5]> satisfies 

r 1 



t»T* &/ (fr.iWT.a, ) 



pr(Tj = 6 I Tj = Of, i / j) oc < 



(b = dj), (6a) 
(n = l,r a ^0,6^aj), (6b) 
w T * tb w T *, aj / (w Tjb w T , aj ) (ri / l,r 2 / 0,6 / a,), (6c) 
^ w T * >h w T * >a ./w T>aj (ri / 1, r 2 = 0, 6 / aj), (6d) 

where r* = (a\, . . . , ctj-i, 6, Oj+i, • • • , &&)• Although r* may not belong to it corresponds 
to a unique partition of and we can use the bijection &k — > 3*1 to recode r* into an element 
of In (l6al)-(l6dl) the event {n = 1, r 2 = 0, 6 7^ aj} is missing since {n = 1, r 2 = 0} implies 
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Conditional simulations of max-stable processes 1 1 

Table 1. Sample path properties of the max-stable models. For the Brown— Resnick 

model, the variogram parameters are set to ensure that the extremal coefficient function 

satisfies 0(115) =1-7 while the correlation function parameters are set to ensure that 

0(100) = 1-5 for the Schlather model. 
Brown-Resnick: -y(h) = (h/\) K Schlather: p{h) = exp{-(/i/A) K } 

71 : Very wiggly 72:Wiggly 73: Smooth pi : Very wiggly p2'. Wiggly p-$: Smooth 
A 25 54 69 208 144 128 

K 0-5 1-0 1-5 0-5 1-0 1-5 

that r* = r, where the equality has to be understood in terms of elements of and this case 
has been already taken into account with (l6ab . 

Once these conditional weights have been computed, the Gibbs sampler proceeds by updat- 
ing e ach element of r successively. We use a random scan implementation of the Gibbs sam- 



pler ( Liu et all 



19951) . More precisely, one iteration of the random scan Gibbs sampler selects 
an element of r at random according to a given distribution, say p = (pi, ■ ■ ■ ,Pk), and then 
updates this element. Throughout this paper we will use the uniform random scan Gibbs sam- 
pler for which the selection distribution is assumed to be a discrete uniform distribution, i.e., 

p = (k-\...,k- 1 ). 



4. Simulation Study 

In this section we check if our algorithm is able to produce realistic conditional simulations 
of Brown-Resnick and Schlather processes. For each model, we consider three different sample 
path properties, as summarized in TableQ] These configurations were chosen such that the spatial 
dependence structures are similar to our applications in Section [5] 
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Fig. 1. Pointwise sample quantiles estimated from 1000 
conditional simulations of max-stable processes with stan- 
dard Gumbel margins with k — 5, 10, 15 conditioning lo- 
cations. The top row shows results for the Brown-Resnick 
models with semi variograms 73 , 72 , 71 , from left to right. 
The bottom row shows results for the Schlather models 
with correlation functions pz,p2,pi, from left to right. The 
solid black lines shows the pointwise 0-025, 0-5, 0-975 
sample quantiles and the dashed grey lines that of a stan- 
dard Gumbel distribution. The squares show the condi- 
tional points {(xi, Zi)}i=i,....fc. The solid grey lines cor- 
respond the simulated paths used to get the conditioning 
events. 

In order to check if our sampling procedure is accurate and given a single conditional event 
{Z(x) = z} for each configuration, we generated 1000 conditional realizations with standard 
Gumbel margins. Figure Q] shows the pointwise sample quantiles obtained from these 1000 sim- 
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Conditional simulations of max-stable processes 13 
Table 2. Computational timings for conditional simulations of max-stable processes on a 

50 x 50 grid defined on the square [0, 100 x 2 1 / 2 ] 2 for a varying number k of conditioning 

locations uniformly distributed over the region. The timings, in seconds, are mean values over 

100 conditional simulations; standard deviations are reported in brackets. 

Schlatter: p(h) = exp {-(>/208) ' 50 } 



Brown-Resnick: y(h) = (h/25) ' 5 
Step 1 Step 2 Step 3 Overall 

fc = 5 0-21 (0-01) 49(11) 1-4(0-1) 50(11) 
fc = 10 8 (2) 76(18) 1-4 (0-1) 85 (19) 

fc = 25 95 (38) 117 (30) 1-4(0-1) 214(61) 
fc = 50 583 (313) 348 (391) 1-5 (0-1) 931 (535) 
Conditional simulations with k — 5 do not use a Gibbs sampler. 



Step 1 Step 2 Step 3 Overall 

1-40(0-02) 1-9(0-7) 0-9(0-3) 4-2(0-8) 

12(4) 2-4(0-8) 1-0(0-3) 15(4) 

86(42) 4(1) 1-0(0-3) 90(43) 

367 (222) 62 (113) 1-0(0-3) 430(262) 



ulated paths and compares them to unit Gumbel quantiles. As expected the conditional sample 
paths inherit the regularity driven by the shape parameter k and there is less variability in re- 
gions close to conditioning locati ons. Since the considered Brown-Resnick processes are ergodic 



( Ka bluchko and Schla ther , 



2010), for regions far away from any conditioning location the sam- 
ple quantiles converges to that of a standard Gumbel distribution indicating that the conditional 
event has no influence. This is not the case for the non-ergodic Schlather processes. Most of the 
time the sample paths used to get the conditional events belong to the 95% pointwise confidence 
intervals, corroborating that our sampling procedure seems to be accurate. 

Table [2] gives computational timings for conditional simulations of max-stable processes on 
a 50 x 50 grid with a varying number of conditioning locations. Due to the combinatorial com- 
plexity of the partition set the timings increase rapidly with respect to the number of con- 
ditioning points k. It is however reassuring that the algorithm is tractable when k 6 {1, . . . , 50}; 
hence covering many practical situations and applications. 
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Fig. 2. Left: Map of Switzerland showing the stations of 
the 24 rainfall gauges used for the analysis, with an insert 
showing the altitude. The station marked with a triangle 
corresponds to Zurich. Middle: Summer maximum rain- 
fall values for 1962-2008 at Zurich. Right: Comparison 
between the pairwise extremal coefficient estimates for the 
51 original weather stations and the extremal coefficient 
function derived from a fitted Brown-Resnick processes 
having semi variogram j(h) = (h/X) K . The grey points 
are pairwise estimates; the black ones are binned estimates 
and the curve is the fitted extremal coefficient function. 
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5. Application 



5 1. Extreme precipitations around Zurich 



In this section we obtain conditional simu lations of extreme pr ecipitation fields. The data 



considered here were previously analyzed by iDavison et al 



(|201 lh who showed that Brown- 



Resnick processes were one of the most competitive models among various statistical models for 
spatial extremes. 

The data are summer maximum rainfall for the years 1962-2008 at 51 weather stations in the 
Plateau region of Switzerland, provided by the national meteorological service, MeteoSuisse. To 
ensure strong dependence between the conditioning locations, we consider as conditional loca- 
tions the 24 weather stations that are at most 30km apart from Zurich and set as the conditional 
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673 Table 3. Distribution of the partition size for the rainfall data esti- 

614 mated from a simulated Markov chain of length 15000 

Partition size 1 2 3 4 5 6 7-24 

Empirical probabilities (%) 66-2 28 4-8 0-5 0-2 0-2 <005 

676 
677 
678 
679 
680 
681 
682 
683 
684 
685 
686 
687 



689 
690 
691 
692 
693 
694 
695 
696 
697 
698 
699 



values the rainfall amounts recorded in the year 2000, the year of the largest precipitation event 
ever recorded in Zurich between 1962-2008, see Figure |2 The largest and smallest distances 
between the conditioning locations are around 55km and just over 4km respectively. 

A Brown-Resnick process having semi variogram ^y(h ) = (h/\) K has to be fitted and 



the maximum pairwise likelihood estimator introduced by IPadoan et all (120 10T) was used to 



simultaneously fit the marginal parameters and the spatial dependence parameters A and 



k. In accordance with 



Davison et al. 



(120111) . the marginal parameters were described by 
7](x) = f3 0jV + /3i jJ? lon(x) + /3 2 ,r;lat(x), a(x) = /3 ,<t + /3i I<T lon(x) + /3 2 , CT lat(x), = /3 ,£, 
where r/(x), <r(x), are the location, scale and shape parameters of the generalized extreme 
value distribution and lon(x), lat(x) the longitude and latitude of the stations at which the data 
are observed. The maximum pairwise likelihood estimates for A and k are 38 (14) and 0-69 
(0-07) and give a practical extremal range, i.e., the distance h + such that 9{h + ) = 1-7, of around 
1 15km, see the right panel of Figure |2] 

Table [3] shows the distribution of the partition size estimated from a Markov chain of length 
15000. Around 65% of the time the summer maxima observed at the 24 conditioning locations 
were a consequence of a single extremal function, i.e., only one storm event, and around 30% of 
the time a consequence of two different storms. Since the simulated Markov chain keeps a trace 
of all the simulated partitions, we looked at the partitions of size two and saw that around 65% of 
the time, at least one of the four up-north conditioning locations was impacted by one extremal 
function while the remaining 20 locations were always influenced by another one. 



16 



C. DOMBRY, F. EYI-MlNKO AND M. RlBATET 



721 
722 
723 
724 
725 
726 
727 
728 
729 
730 
731 
732 
733 
734 
735 
736 
737 
738 
739 
740 
741 
742 
743 
744 
745 
746 
747 








} □ 





1.08 
1.06 
1.04 

- 1.02 

- 1.00 



Fig. 3. From left to right, maps on a 50 x 50 grid of 
the pointwise 0-025, 0-5 and 0-975 sample quantiles for 
rainfall (mm) obtained from 10000 conditional simula- 
tions of Brown-Resnick processes having semi variogram 
j(h) — (ft/38) ' 69 . The right most panel plots the ratio of 
the pointwise confidence intervals with and without tak- 
ing into account the parameter estimate uncertainties. The 
squares show the conditional locations. 



Figure[3]plots the pointwise 0-025, 0-5 and 0-975 sample quantiles obtained from 10000 condi- 
tional simulations of our fitted Brown-Resnick process. The conditional median provides a point 
estimate for the rainfall at an ungauged location and the 0-025 and 0-975 conditional quantiles a 
95% pointwise confidence interval. As indicated by our simulation study, see Figure [T] the shape 
parameter k has a major impact on the regularity of paths and on the width of the confidence 
interval. The value k 0-69 corresponds to very wiggly sample paths and wider confidence in- 
tervals. To assess the impact of parameter uncertainties on conditional simulations, the ratio of 
the width of the confidence intervals with or without parameter uncertainty is shown in the right 
panel of Figure [3] The uncertainties were taken into account by sampling from the asymptotic 
distribution of the maximum composite likelihood estimator and draw one conditional simulation 
for each realization. These ratios show no clear spatial pattern and the width of the confidence 
interval is increased by an amount of at most 10%. 
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Fig. 4. Left: Topographical map of Switzerland showing 
the sites and altitudes in metres above sea level of 16 
weather stations for which annual maxima temperature 
data are available. Middle: Times series of the daily max- 
ima temperatures at the 16 weather stations for year 2003. 
The 'o', ' + ' and 'x' symbols indicate the annual maxima 
that occurred the 11th, 12th and 13th of August respec- 
tively. Right: Comparison between the fitted extremal co- 
efficient function from a Schlather process (solid red line) 
and the pairwise extremal coefficient estimates (gray cir- 
cles). The black circles denote binned estimates with 16 
bins. 

5-2. Extreme temperatures in Switzerland 
In this section we apply our results to get conditional simula tions of extreme temperature 



fields. The data considered here were previously analyzed by 



Davison and Gholam-Rezaee 



(1201 1[) and consist in annual maximum temperatures recorded at 16 sites in Switzerland dur- 



ing the period 1961-200 5, seeFigu re[H 



Following the work of 



Davison and Gholam-Rezaed (12011 



), we fit a Schlather process with an 



isotropic powered exponential correlation function and trend surfaces r](x) = (3o tV + /3i ir) alt(x), 
a(x) = /3o,o-> C( x ) = Po,£ + /3i,£alt(x), where alt(x) denotes the altitude above mean sea level 
in kilometres and {r]{x), a(x), are the location, scale and shape parameters of the gener- 

alized extreme value distribution at location x. The spatial dependence parameter estimates are 
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Table 4. Distribution of the partition size for the temperature 

data estimated from a Markov chain of length 10000 
Partition size 1 2 3 4 5-16 

Empirical probabilities (%) 2-47 21-55 64-63 10-74 0-61 



821 

A = 260 (149) and k = 0.52 (0.12) and the corresponding fitted extremal coefficient function, 

822 

similar to some extent to our test case p% in Section HI is shown in the right panel of Figure [H 

823 

In year 2003, western Europe was hit by a severe heat wave believed to be the hottest one ever 

824 

recorded since at most 1540 ("2003 European heat wave", Wikipedia: The Free Encyclopedia). 

825 

Switzerland was largely impacted by this severe extreme event since the nation wide record tem- 

826 

perature of 41-5°C was recorded that year in Grono, Graubunden, near Lugano. Consequently 

827 

for our analysis we use as conditional event the maxima temperatures observed in summer 2003, 

828 

see FigurelU Based on the fitted Schlather model, we simulate a Markov chain of effective length 

829 

lOOOOwith a burn-in period of length 500 and a thinning lag of 100 iterations. The distribution 

830 

of the partition size estimated from these Markov chains is shown in Table [4j We can see that 

831 

around 90% of the time the conditional realizations were a consequence of at most three ex- 

832 

tremal functions. Since our original observations were not summer maxima but maximum daily 

833 

values, a close inspection of the times series in year 2003 reveals that the hottest temperatures 

834 

occurred between the 11th and 13th of August, see Figure |4] and, to some extent, corroborates 

835 

the distribution of Table [4] 

836 

Figure [5] shows the 0-025, 0-5 and 0-975 pointwise sample quantiles obtained from 10000 

837 

conditional simulations on a 64 x 64 grid. As expected, we can see that the largest temperatures 

838 

occurred in the plateau region of Switzerland while temperatures were appreciably cooler in 

839 

the Alps. The right panel of Figure [5] shows the difference between the pointwise conditional 

840 

medians and the pointwise unconditional medians estimated from the fitted trend surfaces. The 

841 
842 
843 
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Fig. 5. From left to right: Maps on a 64 x 64 grid of the 
pointwise 025, 0-5 and 0-975 sample quantiles for tem- 
perature (°C) obtained from 10000 conditional simula- 
tions of the fitted Schlather process. The squares show the 
conditional locations and the associated conditional values. 
The right most panel shows temperature anomalies, i.e., the 
difference between the pointwise conditional and uncondi- 
tional medians. 

differences range between 2-5° C and 4-75° C and the largest deviations occur in the plateau 
region of Switzerland. 
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Appendix 

The Brown-Resnick model 
For all x € X k and Borel set Act' 

A X (A)= / pr[(exp{W(x)-j{x)}eA}C 2 d(= / / l{c^ P { y - 1 ( x )}eA}f x (y)dyC 2 dC, 

JO JO JR k 

where f x denotes the density of the random vector W(x), i.e., a centered Gaussian random vector with 
covariance matrix Yj x and variance 27(3;). The change of variables z — (exp{y — 7(2;)} and r = log£ 
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913 yields 



914 
915 
916 



936 
937 
938 
939 



/CO n n 
/ f x {\ogz-r + 1 (x)}\{z- 1 dze- r dr= / \ x (z)d 
-go J A J A 



with 



/>00 

917 ^) = H^ 1 / x {log^-r + 7(a;)}e- r dr. 

„-_i J -co 



Since 



918 
919 
920 

921 with 



f x {logz - r + 7 (z)}e- r - {2n)- k ' 2 \i: x \- 1 / 2 cxp , 



922 
923 
924 
925 
926 
927 
928 

929 and since log(u, z) — J m ^ log u + J m ^ log z, it is not difficult to show that 

930 
931 
932 
933 

934 The Schlather model 

935 For all x e X k and Borel set AcR k 



P(r) = rHl^-Hu - 2r [l^E^log* + 7 (x)} - l] + {logz + 7 ( a ;)} T £- 1 {log z + 7 (x)}, 

standard computations for Gaussian integrals give 

X x (z) = C x exp ^-ilogz T Q x logz + L^logz^ fl^ 1 - 

The conditional intensity function is 

C ( 1 1 1 m 



C ( 1 1 m 

A 5 |x, 2 (w) = -^■exp|--(logu-/i s | a; ^) T S^(logu-/i iS | a; ,j| 11" * 

Finally, the relation C( s , x )/C x = (27r) _m / 2 |£ s | x | -1 / 2 is a simple consequence of the normalization 

fK\x,zi u ) du = L 



/•'DC 1 pOO p 
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where f x denotes the density of the random vector s(x), i.e., a centered Gaussian random vector with 



covariance matrix E x . The change of variable z — V^Cu gives 

- {2^)- k \Y lx \- 1/2 I -^ =r -E[X k ~ 1 ]dz, X ~ Weibull ( 

= (27T)- fe |S x |- 1 / 2 J 



Air 



Z^Yj x z 



4, ^-«V*+i 1<l; 



2tt 

K(z)dz, 

where A,(z) = tt"^" 1 )/ 2 |S,|- 1 / 2 a :c (z)-( fe + 1 )/ 2 r {(fc + l)/2} and a x (*) = ^E" 1 ^ 
For all u € M" 1 the conditional intensity function is 



K\x,z{ u ) = n 



V2 l S (^)l V2 \ a (s,x){u,z) 



-(m+fe+l)/2 



-ro/2_ 



'm+fc+1 ' 



sk,zV ^ " |E X |-V2 i a x (s) 

We start by focusing on the ratio a^ s ^(u, z)/a x (z). Since 

(E s — £ s:x £ x 1 Yi x:s ) 1 — (E s — E s;x E x 1 E a;:s ) 1 S s;x S x 1 

T x E x:s (E s — £ s;x E x E x;s ) 1 E x 1 + E x E x:s (E s — £ s;x E x £ x;s ) 1 Yi s - x Yi x 1 
straightforward algebra shows that 





-1 i- 


^s:x 









a x (^) fc + 1 

We now try to simplify the ratio |E( S;X J/|E X |. Using the fact that 



^s:x 




Id m Xl s::c 






^X'.S ^x 











combined with some more algebra yields 

— p=— j — — |Zj s Zj s - x Zj ^x:s\ 



k+ 1 

a K (z) 



ISI. 



Using the two previous results it is easily found that 

™ ~ , ^ -(m+fc+l)/2 . . . , , . 

(ti-zifs-Hu-A*)! r(H^±i) 



A fl |^( U )-7r- m /2(fe + l)-W2|f;|-i/2 1 + 



fe + 1 



r(¥) ' 
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which corresponds to the density of a multivariate Student distribution with the expected parameters. 
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